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In a recent experiment, signatures of Majorana fermion (MF) were found in the vortex core threading a 
heterostructure composed of n layers of topological insulator (TI) deposited on a bulk s-wave 
superconductor. Here we provide strong theoretical support to the experiment. First, we demonstrate that 
MF modes appear on both top and bottom layers of TI, and are well separated for n ~ 6. The top MF becomes 
more extended with increasing n, in agrement with the experiment. Second, we show both analytically and 
numerically that right at the vortex core the MF mode is always accompanied by another low energy bound 
state, leading to a zero-bias peak plus a side peak in the local density of states (LDOS) therein. However, a 
local scalar impurity at the core can wipe out the accompanying side-peak state while leaving the zero-energy 
MF mode intact. Consequently the LDOS becomes symmetric about the fermi level, and the peak does not 
branch near the vortex core, in agreement with the experiment. Finally but unfortunately, while the MF is 
extremely stable against a single local impurity, the stability in terms of the critical impurity strength is 
reduced drastically for a moderate concentration (e.g., 10%) of impurities. 

A Majorana fermion (MF) is an anti-particle of itself, originally proposed as a candidate for neutrino 1 . 
Recently it was proposed that MF may be generated in some particular condensed matter systems. For 
example, a zero energy MF mode can appear in the vortex core of the s-wave superconductor/bulk 
topological insulator heterostructure 2 . As a variant, MF may also appear at the open ends of a semiconducting 
nano-ribbon with strong spin-orbital coupling, in proximity to an s-wave superconductor and under a high 
magnetic field 3 . Apart from such hetero-structures, MF may appear in vortices of pure systems, such as the 
superfluid 3 He-B 4 and the potential + ip y -wzve superconductor Sr 2 Ru0 4 5 . The renewed extensive interests in 
search of MFs are mainly driven by the potential application of MFs. Since the MFs obey non-Abelian statistics, 
they can be utilized to perform fault-tolerant topological quantum computing 6 . In addition, the MF mode also 
leads a plethora of new intriguing effects, such as fractional Josephson effect 7 ' 8 , resonant Andreev reflection 910 , 
and resonant crossed Andreev reflection 1112 . 

Recently important experimental progresses have been made for the topological insulator/s-wave supercon- 
ductor heterostructure 1314 . In the earlier theoretical proposal 2 , a thin layer of superconductor (SC) is to be 
deposited on a bulk topological insulator (TI). However, experimentally it is difficult to prepare a pristine bulk 
TI (such as Bi 2 Se 3 ), while it is much easier to prepare a bulk SC (such as Ne 2 Se 3 ). Therefore, in Ref. 13 a few 
(quintuple) layers of TI was deposited on the bulk SC substrate, resulting in a TI/SC configuration. Signatures of 
both superconductivity and the surface states of TI (away from the Fermi level) were observed 13 . More recently, 
the structure is further improved by replacing Bi 2 Se 3 with Bi 2 Te 3 14 . Scanning tunneling microscopy measure- 
ments (STM) of the local density of states (LDOS) on the top layer revealed that the zero-bias peak in the LDOS is 
perfectly symmetric and does not branch in energy up to a finite distance away from the vortex core, if the number 
of TI layers n > 4. This non-branching behavior is in contrast to the case of n < 4, in which the branching begins 
right away from the vortex core as in usual s-wave superconductors. The non-branching behavior was argued as 
the signature of MF. The experiment appears very promising. However, in order for the heterostructure to be 
applied, a few important issues need to be addressed. First, since the zero-bias-peak appears for any case of n in the 
experiment, it is necessary to identify/falsify the zero energy bound state as a MF, in addition to the identification 
based on the non-branching behavior. Second, it is necessary to determine and understand the optimal number of 
TI layers to realize MF. Finally, in a more general respect, there are inevitably imperfections or impurities in real 
samples. For applications it is important to understand how impurities would destroy a MF mode. We should 
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Figure 1 | (a) Schematic plot of a vortex threading the IX IX n TI above the SC substrate (gray box), (b) The assumed proximity function j{z). 



mention that according to the experiment, the Fermi level is away 
from the Dirac point of the TI surface states (in fact close but below 
the conduction band). The situation is very different to the SC/TI/SC 
configuration with the Fermi level right at the Dirac point addressed 
elsewhere 15 " 17 . 

Here we address the above issues theoretically. First, we dem- 
onstrate that MF modes appear on both top and bottom layers of 
TI, and are well separated for a moderate n. The top MF becomes 
more extended with increasing n, in agrement with the experiment. 
Second, we show both analytically and numerically that right at the 
vortex core the MF mode is always accompanied by another low 
energy bound state, leading to a zero-bias peak plus a side peak in 
the local density of states (LDOS) therein. However, a local scalar 
impurity at the core can wipe out the accompanying side-peak state 
while leaving the zero-energy MF mode intact. Consequently the 
LDOS becomes symmetric at the vortex core, and the peak does 
not branch near the vortex core, in agreement with the experiment. 
Finally but unfortunately, while the MF is extremely stable against a 
single local impurity, the stability in terms of the critical impurity 
strength is reduced drastically for a moderate concentration of 
impurities. 

The 3d heterostructure. The in-plane structure of the experimental 
TI is a triangle lattice, and each unit cell contains a quintuple layer in 
the out-of-plane direction. For convenience but without loss of 
generality, we model the 3D TI by a simple cubic lattice, with the 
hamiltonian 

Ho= 1 Y [^ t tb(r 0 -/r b ).A 1 . +b +h.c.] 



2 1> 



(i) 



fa- 



Here i// r is a four- spinor annihilation field at site r = (x,y,z),b = x,y,z 
is a bond vector connecting nearest-neighbor sites, and t h is the 
hopping parameter along bond b. On the other hand, T 0 = rj 3 0 cr 0 
and r b = rji (x) a b are four Dirac matrices. The Pauli matrices r\ (a) 
act on orbital (spin) basis. The parameters m controls the topology of 
the insulator, and finally \i is the chemical potential. We take z as the 
direction normal to the quintuple layer in Bi 2 Se 3 or Bi 2 Te 3 , and we 
model the anisotropy by taking (t^tyJz) = (t\\,t\\,t±). 

We assume that the SC substrate contacting the bottom layer leads 
to an effective s-wave pairing hamiltonian 

ff p =£#A(r)fo># + h.c. (2) 

r 

Under an applied magnetic field along z, vortices appear and thread 
the heterostructure. Experimentally the vortices are dilute for a mag- 



netic field of the order of Tesla. Thus we consider an isolated vortex 
line in an L X L X n cubic lattice, as schematically shown in Fig. 1(a). 
The pairing potential is written as, 



A(r) = A 0 taiih(r l /{)e'ty(z) 



(3) 



where A 0 is the pairing amplitude in SC, rj | = (x, y) is the planar 
position vector (the origin of which is set at the center of the L X L 
lattice) with azimuthal angle <j> , £ is the coherence length of the SC 
substrate, and finally f(z) is a proximity function. The concrete form 
oif(z) is unknown. However it is reasonable to take /(h) = A n /A 0 , 
where A n is the pairing gap measured at the top layer of a hetero- 
structure with n layers of TI. In qualitative accordance to the experi- 
ment 13,14 , we take such a function as plot in Fig. 1(b). We assume £ = 
£o/f(z), where £ 0 is the coherence length in the SC substrate, to reflect 
the fact that the effective pairing gap varies along z. In the London 
limit the magnetic field leads to a vector potential a(r) = n(—y, x, 0)/ 
(2L 2 ) in dimensionless units. This amounts to a unit of flux (j) 0 = he/ 
2e threading the planar area. The vector potential couples to elec- 
trons minimally via 

Jb-a(r+b/2) 



t h ^t h e l 



(4) 



in the first line of Eq.(l). 

The complete hamiltonian is H = H 0 + H p . Formally we rewrite it 
in the Nambu space as 



(5) 



with the spinor ^ = (fa^ r ), and the corresponding single-particle 
hamiltonian h(r, r') (each being an 8 X 8 matrix acting on spin, 
orbital and Nambu bases). It is easy to show that h satisfies the 
anti- commutation relation {tiK, h} = 0, where x x is the first Pauli 
matrix acting in the Nambu space and K is the complex conjugation 
operator. By this antisymmetry, iih\ x ¥) = it is guaranteed that 
h^K^) = -EtiK|¥>. Thus |¥) and t^Y) form a pair of (ortho- 
gonal) states with opposite energies. Because of the doubling of 
degrees of freedom in the Nambu space, only one of such two states 
are independent. If E = 0, the two states become degenerate, and any 
linear combination of them is also an eigenstate. In a particular 
combination resulting in a pair of new states l*?^), it is possible to 
enforce t^I^) = and t 1 K\ x ¥ 2 ) = |¥ 2 > (up to trivial phases), in 
addition to the orthogonality condition (^l^) = 0- The quasi- 
particle field operators associated with these states are self-hermitian, 
defining the MF modes. They would be stable against perturbations 
provided they are spatially bounded and not overlapping each other. 
Such a MF mode was argued to appear in the vortex core in a hetero- 
structure with a layer of s-wave SC on a bulk TP. The situation we are 
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Figure 2 | The wavefunction profile p 0 (r) for (a) n — 2, (b) 12 = 4, (c) 12 = 

6, and (d) a = 8. In each panel, a corner is cutoff to reveal the interior of the 
TL 



addressing is upside down in structure: n layers of TI on top of a bulk 
s-wave SC. For sufficiently large n, helical states appear on the sur- 
faces of the otherwise isolated TI. By coupling to SC, there should be a 
MF mode near the interface (the bottom layer of TI) according to Ref. 
2. For the same reason another MF mode should appear on the top 
layer if the proximity induced pairing gap persists over there. For an 
intermediate n, the postulated MF's near the bottom and top layers 
may couple to each other, lifting the energy degeneracy. 

To proceed we investigate numerically the evolution of the lowest- 
lying state I^Fq) (with energy E 0 closest to zero) versus n to gain 
insights into how MF states would emerge out of |^ 0 ) ana 
TxKl^Po). We define the probability density function 



to reveal the detailed spatial profile of the emerging MF wavefunc- 
tion^) in the heterostructure. Here the trace is taken over spin, 
orbital and Nambu degrees of freedom. We notice thatp 0 ( r ) is ident- 
ical for both |*F 0 ) and Ti-fCj^o), although they are necessarily ortho- 
gonal if E 0 7^ 0. In the calculations we fix the parameters (t\\, t± , m, fi, 
A 0 , i 0 ) = (1, 0.5, 0.3, 0.23, 0.1, 10) without lost of generality. The 
smaller hopping t ± takes qualitative account of the anisotropy in real 
samples. Fig. 2 shows the resulting p 0 (r) (color scale) for n = 2, 4, 6, 8 
and L = 150 = 15cj 0 . We see that the state is bounded to the vortex 
line in all cases. It is fully connected across the layers for n = 2 in 
Fig. 2(a). In Fig. 2(b), p 0 (r) is mainly concentrated near the top and 
bottom layers, implying that this is not the usual Andreev bound 
state, but rather two MF states slightly reconstructed due to the 
coupling in the interior of the sample. Indeed, a separate calculation 
shows that in an isolated n > 4 sample, gapless helical surface states 
begin to show up. Such states coupled to the pairing potential would 
lead to MF states. For larger n in Fig. 2(c) and (d), the function p 0 (r) 
becomes more disconnected along z due to weaker coupling between 
the top and bottom layers, and at the same time more extended on 
the top layer where the effective pairing gap decreases for increasing 
n. Since an extended MF mode would be fragile, the above results 
suggest that there is an optimal value of n, say n opti to look for MF 
mode on the top layer of the present heterostructure. Since the energy 
E 0 is nonzero in our finite-sized sample, we check the scaling of E 0 
versus 1/L in Fig. 3. Clearly E 0 is finite and asymptotically independ- 
ent of L in the limit of large L for n = 2, while E 0 —> 0 with increasing 
L for n = 6. The case for n = 4 appears marginal due to the coupling 
between top and bottom layers. Thus we expect a genuine MF mode 
on the top layer for n > n opt = 6 in the thermodynamic limit of L, 
supporting the experimental setup in search of MF. Of course the 
above n opt should only be taken qualitatively, given the difficulty to 
extract the exact parameters in the experiment. 

Before closing this section we remark that since the top and bot- 
tom MF's are bound to the vortex line, braiding two vortices leads to 
braiding of the MF pairs. The statistics during such a braiding is 
identical to that for conventional fermions. Thus in the present het- 
erostructure it might be difficult to achieve non-Abelian statistics 
that requires braiding individual MF's. 
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1/L 

Figure 3 | Finite size scaling of E 0 versus 1/L for various n. Symbols are the data and lines are guides to the eye. 
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Effective 2d models. In the experiment the LDOS is measured by 
STM. Since the contribution mainly comes from the top layer, it is 
beneficial to investigate the top layer in more details. To reduce the 
computation load, we resort to an effective 2D model describing the 
low-lying quasi-particle states. We first gain understanding by 
analytical analysis of the a continuum model. We then perform 
numerical simulations in a lattice model. We also study the 
robustness of the MF state against impurities, and compare the 
LDOS to the unexpectedly symmetric zero bias peak observed in 
the experiment. 

Analytical treatment in the continuum limit: In the low energy 
window, the effective hamiltonian in the Nambu space can be written 
as, 



H = J [( - tV-a -h)t z - A^ x + A 2 T y ] 



(?) 



Here we used a special convention for the annihilation field 
*¥= (tf/pil/^, — ^J,i^|^ > and the group velocity v D of the helical 
states is set to unity. Henceforth we use r to denote the 2D position. 
We performed a spin rotation to write V * a instead of d x o y — d y o x * 
Ai (A 2 ) is the real (imaginary) part of the s-wave pairing potential 

A(r)=Ai(r) + iA 2 (r)=A 0 tanh(r/{)<* (8) 

where A 0 (r) is the pairing scale, £ the coherence length and </> the 
azimuthal angle of r. Finally we did not include the vector potential a 
since for an isolated vortex in an infinite area and in the London limit 
the effect of a is negligible as far as the low energy bound states are 
concerned. In the normal state, the hamiltonian leads to the band 
dispersion E k = ± k — /n, where the sign depends on the chirality ha. 
Because of the spin-momentum locking the two bands are spin non- 
degenerate. Therefore the intra-band pairing must be odd, and turns 
out to hep + /p'-like in the present case. This implies immediately a 
MF in the vortex core, just as in a p + ip' -superconductor 18 . The 
inter-band pairing is irrelevant for A 0 <C|ju|. 

For our purpose, we need the wavefunctions of the low lying states. 
Following Volovik 18 , this is most easily obtained in the quasi- classical 
approximation. By rotation symmetry the states can by classified by 
the angular momentum /. The energy up to the first order of / is E = 
— l(D 0 , with 



OJQ 



= Jdsf 0 2 (s) 



A 0 tanh(r/g) 
k F r 



(9) 



where k F is the Fermi momentum and/ 0 is an envelop function 
1 



/o(s)~-exp 



- f A 0 tanh(r7£) 
J o 



cosd'ds' 



(10) 



with a normalization constant C and s = r cos (j). The corresponding 
wavefunction is 



%(r)~2ni%(r) 



iJ l+1 (k P r)e>« + ^ 

h.^kpryv-w 

\ ihikpry* ) 



(11) 



where // is the Bessel function. The energy levels obtained here are 
essentially similar to the Caroli-de Gennes-Matricon states in usual 
s-wave superconductors 19 , except that the angular momentum / here 
is an integer rather than a half integer due to the fact that the helical 
state picks up a negative sign upon a rotation by In. 

Of particular interest is the state with / = 0. The general argument 
presented in the previous section implies this zero energy state must 
describe a MF mode. Indeed, the creation field operator associated 
with this state is, up to the envelop function f 0 (r), 



~/o«/4 - iM^ + ihe% +J-ie~%, 



(12) 



which is self-hermitian up to a trivial global phase. In principle there 
should be another MF mode on the outer boundary of an open 
system. But this is an artifact of the 2D model in view of the 3D 
heterostructure (where MF's appear on the top and bottom layers). 
The electronic LDOS from the above states is given by 

p(co,r)x [// 2 (M+// 2 -i(M]/o(r)^(a;-/a;o). (13) 



Since //(0) = 0 for / # 0, the contribution to the spectral weight at r = 
0 comes from / = 0 and / = 1 at energy co = 0 and w = w 0 , 
respectively. Thus there are two LDOS peaks located asymmetrically 
in energy. The energy of the side peak is positive for /n > 0 here, but 
would be negative if fi < 0 instead. The asymmetry in the LDOS is a 
direct consequence of the absence of the particle-hole symmetry in 
the normal state band structure since the Fermi level is away from the 
Dirac point in our case. Similar asymmetries were actually found 
even for the bound states in the vortex core of usual s-wave super- 
conductors as the Fermi level is close to the band edges, both the- 
oretically and experimentally 20 ' 21 . Of course, as the Fermi level 
approaches the Dirac point (/n — > 0), k F decreases and thus the level 
spacing w 0 in our case increases. As k F —> 0, the quasi-classical 
approximation fails but the tendency of co 0 oc l/k F does imply the 
absence of in-gap bound states other than the unique Majorana zero 
mode for k F = 0, as discussed elsewhere 15 " 17 . But this is clearly not the 
case in the experiment, where many bound states are observed 14 . 

As we step away from the vortex core, the other states with higher / 
contributes, and from the behavior of Bessel functions, LDOS peaks 
should appear at the guiding centers of Jv±i(k F r) for co = ± lw 0 . Such 
a peak splitting behavior is similar to that in usual superconductors. 
Rather unexpectedly the experiment reports a symmetric zero -bias 
peak, and it appears that the central peak does not split until r is a few 
atoms away from the vortex core 14 . One possible reason for the 
discrepancy is the working temperature is higher than co 0 that smears 
out the asymmetry due to the side peak. Another possibility is impur- 
ities pinning the vortex may have suppressed the / = 1 state locally, 
leaving the / = 0 Majorana state intact. We shall come back to this 
point later. 

Lattice model: We now consider an effective 2D lattice model for 
the low-lying states on the top layer of the heterostructure. The 
effective hamiltonian is, 



r,b 



(14) 



where b = x,y. We set t h = t = 1 and /n = 0.23. Notice that the o z 
terms are necessary to make sure that there is only one Dirac point in 
the Brillouine zone. Such terms violate time-reversal symmetry. 
Nonetheless, this effect is negligible for low lying states close to the 
Dirac point. In the presence of vortex, we couple the vector potential 
to electrons by minimal coupling as described in the previous section. 
We use a large lattice (of linear dimension L = 400 in this subsection) 
to reduce finite size effect. The pairing hamiltonian H p is of a form 
closely similar to that in Eq.8, except that r is a 2D vector here. From 

now on we switch back to the convention ^ = in 

the Nambu space. The LDOS is given by 



p(r,co)= ^2S(co-E n )A n (r), 



(15) 



where 



SCIENTIFIC REPORTS | 4:6363 | DOI: 1 0.1 038/srep06363 



4 



A„(r)=Tr^± Ii |«P„(r))(«P„(r)|, 



(16) 



is the level-resolved local spectral weight. In practice, the delta func- 
tion is smeared by a proper factor. On the other hand, the dimension 
of the single particle hamiltonian is too large to get all eigenstates. 
Instead, we use the real-space recursion method 22 ' 23 to get p(r, co) 
most reliably. 

Fig. 4(a) shows p(r, w) for r along a (1, l)-cut emanating from the 
vortex core, and (b) shows a few low energy levels with the height 
proportional to A n (x). There is clearly a ZBP and a side peak near the 
vortex core. As discussed previously, the zero energy state must be a 
MF mode. The energy of the side peak is E Y ~ A 0 /k F ^ where k F ^ — 2.3 
in our case. Moreover, the peaks shift to higher energies as r is more 
distant from the core. All of these features are consistent with the 
previous analytical discussions. In particular, there should be a side 
peak accompanying the ZBP. It is however not observed in the 
experiment. Apart from thermal smearing, the discrepancy is likely 
caused by impurities which we consider in the following. 

The impurity hamiltonian can be written as, 



Himp — ^ 



J2 V i (r)<r i T z + V 2 (r)<r > 



TO 



(17) 



where Vo(r) and V/= 1,2,3(1*) are the position dependent impurity 
strengths of scalar and magnetic impurities (of various polariza- 
tions), respectively. Here the single-particle part is denoted as 
V imp (r) for convenience. We observe that {V imp , TiK] = 0. For a 



MF state |Y> (satisfying t^Y) 
to the first order of V ir 



I^F)), the energy correction is, 



' imp) 



AE= (V\V imp \V) = <¥| V^nKl*?) 

= -(V\T 1 KV imp \V) = -{V\V imp KT 1 \V) 
= -{V\V imp \V), 



(18) 



which must be zero by comparing the first and last lines. Therefore 
the MF state is immune to weak impurities of any kind. In this sense 
the MF is both charge and spin neutral. On the other hand, no 
protection mechanism applies for the energy of other finite- energy 
bound states against V imp . Fig. 5 is a similar plot to Fig. 4 but with a 
local scalar impurity right located at the vortex core. The impurity 
strength is V 0 = 15A 0 . (This is a weak impurity since A 0 is a small 
energy scale as compared to the bandwidth.) We observe that the 
previous side peak at energy E x is completely wiped out. The higher 
energy levels remain intact since their wavefunctions vanish at the 
vortex core, thus are decoupled to the local impurity to the leading 
order. The LDOS becomes much more symmetric at the vortex core, 




(b) 
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Figure 5 | The same plot as Fig. 4 except that a local scalar impurity with 
Vo = 15A 0 is located at the vortex core. 

and the peak does not break near the vortex core, in nice agreement 
with the experiment 14 . We should emphasize, however, that the 
LDOS in our case is still slightly asymmetric near the gap edges, as 
shown in Fig. 5b. The asymmetry is an intrinsic feature of the model 
as the Fermi level is away from the Dirac point, as we discussed 
analytically, and is enhanced by the relatively large gap scale in the 
calculation that is necessary to improve the numerical resolution. In 
the experiment, the gap scale is much smaller than the Fermi energy 
(the distance from the Fermi level to the Dirac point), so that the 
energy level spacing of the bound states is smaller. Combined with 
the finite temperature smearing, the experimental LDOS may appear 
more symmetric than what we show in Fig. 5. Therefore the differ- 
ence to the theory is quantitative rather than qualitative. 

We also show the corresponding results for a local magnetic 
impurity in Fig. 6. The impurity reduces the amplitude of the MF 
mode at the core, but the zero energy level remains robust. On the 
other hand, additional in-gap states are induced by the impurity. The 
LDOS is asymmetric in energy due to the side peaks. Experimentally 
neither the reduced ZBP nor the asymmetry of LDOS is observed, 
implying the absence of magnetic impurities therein. 

Finally we examine the effect of a finite concentration of impur- 
ities, an issue also relevant to experiments. Since the impurities too 
afar from the vortex core are irrelevant for the low-lying bound 
states, we put impurities within a radius R = 10 around the core 
randomly, with ensemble averages 

(V,(t)V j (r')) = Wf8(R-r)d t sd v , (19) 

where W t is the strength of type i. For magnetic impurities we define 
W m = ^/y^._ T ~ ~ W? as a measure of the total strength. By mon- 
itoring the presence/absence of the zero energy level, we obtain a 
phase diagram in Fig. 7, showing the critical strength W c of W 0 (W m ) 
for scalar (magnetic) impurities, above which the MF disappears, as a 
function of the impurity concentration c. The critical strength is 
reduced to the order of gap scale for both scalar and magnetic impur- 
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Figure 4 | (a) LDOS along a (1, l)-cut, starting from the vortex core 
(bottom line), (b) Low energy levels with the height proportional to the 
respective spectral weight for each point along the cut. 
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Figure 6 | The same plot as Fig. 4 except that a local magnetic impurity 
with strength V x = 12A 0 is located at the vortex core. 
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Figure 7 | W c versus c for scalar (red) and magnetic (blue) impurities. 

ities. This is in contrast to the stability in the single impurity case. The 
result suggests that high sample quality is necessary. 

Summary 

In conclusion, we provided strong theoretical support to the experi- 
mental result that there exists an optimal number (n opt ) of TI layers 
for MF to appear on the top layer of the TI/SC heterostructure. We 
point out however that the MF modes on the top and bottom TI 
layers are bound to the vortex line and would cause difficulty to 
achieve braiding individual MF's. We show analytically and numer- 
ically that the LDOS near the vortex core has a ZBP and a side peak. 
The absence of the side peak in the experiment is nicely reconciled by 
the effect of a local scalar impurity pinning the vortex core. We also 
examined the effect of a finite concentration of impurities. While the 
MF mode is extremely robust to a single local impurity, the critical 
impurity strength drops drastically to the order of the effective pair- 
ing gap at a moderate impurity concentration. This points toward the 
importance of sample quality. We suggest that the effect of single and 
many impurities be tested by STM (at sufficiently low temperatures) 
intentionally to corroborate the MF mode further and to examine its 
stability, paving the way toward the application of such fascinating 
objects. 
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